# Replication: The unequal effect of economic development on perceived labor market risks and welfare
# Tabea Palmtag, University of Zurich
# 24.07.2023

rm(list=ls(all=T))
options(stringsAsFactors = FALSE)

p_needed <- c("brms", "rstan", "dplyr","bayestestR", "stargazer", "lfe")
lapply(p_needed, require, character.only = TRUE)


source("fun_tables.R")

### STAN MODELS ###
rstan_options(auto_write = TRUE)
options(mc.cores = 4)

dat <- readRDS("cleandata_replication.Rds")

## Main Models----
# Labor market risk 
econ30 <- brm(probecon | weights(combinwt) ~ educ * light_30_log + age + female + 
                urban + unempl + ethfair + media +
                conflict + ea_service + ea_infrastructure + as.factor(cyear) +
                (1 | EA), data = dat, family=binomial("logit"),
              chains=4, seed=123, iter=2000, warmup=1000)
#saveRDS(econ30, "fit/econ30.Rds")
mytable(econ30, livingcond = F)

# Living Conditions
liv30 <- brm(as.ordered(plivcond) | weights(combinwt) ~ educ * light_30_log + age + female + 
               urban + unempl + ethfair + media +
               conflict + ea_service + ea_infrastructure + as.factor(cyear) +
               (1 | EA), data = dat, family=cumulative("logit"),
             chains=4, seed=123, iter=2000, warmup=1000)
#saveRDS(liv30, "fit/liv30.Rds")
mytable(liv30, livingcond = T)

## Non-interacted models----
econ30ni <- brm(probecon | weights(combinwt) ~ educ + light_30_log + age + female + 
                  urban + unempl + ethfair + media + 
                  conflict + ea_service + ea_infrastructure + as.factor(cyear) +
                  (1 | EA), data = dat, family=binomial("logit"),
                chains=4, seed=123, iter=2000, warmup=1000)
#saveRDS(econ30ni, "fit/econ30_ni.Rds")
mytable(econ30ni, livingcond = F)

liv30ni <- brm(as.ordered(plivcond) | weights(combinwt) ~ educ + light_30_log + age + female + 
                 urban + unempl + ethfair + media + 
                 conflict + ea_service + ea_infrastructure + as.factor(cyear) +
                 (1 | EA), data = dat, family=cumulative("logit"),
               chains=4, seed=123, iter=2000, warmup=1000)
#saveRDS(liv30ni, "fit/liv30_ni.Rds")
mytable(liv30ni, livingcond = T)

## Labor market risks alternative operationalizations----
econ30_all <- brm(prob_all | weights(combinwt) ~ educ * light_30_log + age + female + 
                urban + unempl + ethfair + media +
                conflict + ea_service + ea_infrastructure + as.factor(cyear) +
                (1 | EA), data = dat, family=binomial("logit"),
              chains=4, seed=123, iter=2000, warmup=1000)
#saveRDS(econ30_all, "fit/econ30_all.Rds")
mytable(econ30_all, livingcond = F)

## Comparison to national development levels // 30km radius----
econ30_diff <- brm(probecon | weights(combinwt) ~ educ * light30_diff + age + female + 
                urban + unempl + ethfair + media +
                conflict + ea_service + ea_infrastructure + as.factor(cyear) +
                (1 | EA), data = dat, family=binomial("logit"),
              chains=4, seed=123, iter=2000, warmup=1000)
#saveRDS(econ30_diff, "fit/econdiff.Rds")
mytable(econ30_diff, livingcond = F)

liv30_diff <- brm(as.ordered(plivcond) | weights(combinwt) ~ educ * light30_diff + age + female + 
               urban + unempl + ethfair + media + 
               conflict + ea_service + ea_infrastructure + as.factor(cyear) +
               (1 | EA), data = dat, family=cumulative("logit"),
             chains=4, seed=123, iter=2000, warmup=1000)
#saveRDS(liv30_diff, "fit/livdiff.Rds")
mytable(liv30_diff, livingcond = T)

## Changing local development----
# only including observations with local growth measures: 3 year change in local development (night light)
change <- dat %>%
  filter(year > 2002)
econ30_ch <- brm(probecon | weights(combinwt) ~ educ * change3_30_log + age + female + 
                urban + unempl + ethfair + media +
                conflict + ea_service + ea_infrastructure + as.factor(cyear) +
                (1 | EA), data = change, family=binomial("logit"),
              chains=4, seed=123, iter=2000, warmup=1000)
#saveRDS(econ30_ch, "fit/econ3ch.Rds")
mytable(econ30_ch, livingcond = F)

liv30_ch <- brm(as.ordered(plivcond) | weights(combinwt) ~ educ * change3_30_log + age + female + 
               urban + unempl + ethfair + media + 
               conflict + ea_service + ea_infrastructure + as.factor(cyear) +
               (1 | EA), data = change, family=cumulative("logit"),
             chains=4, seed=123, iter=2000, warmup=1000)
#saveRDS(liv30_ch, "fit/liv3ch.Rds")
mytable(liv30_ch, livingcond = T)

## Different sources of nightlight----
# only including observations without onshore petroleum fields
petro <- dat %>%
  filter(petro == 0)
econ30_petro <- brm(probecon | weights(combinwt) ~ educ * light_30_log + age + female + 
                      urban + unempl + ethfair + media + 
                      conflict + ea_service + ea_infrastructure + as.factor(cyear) +
                      (1 | EA), data = petro, family=binomial("logit"),
                    chains=4, seed=123, iter=2000, warmup=1000)
#saveRDS(econ30_petro, "fit/econ30_petro.Rds")
mytable(econ30_petro, livingcond = F)

liv30_petro <- brm(as.ordered(plivcond) | weights(combinwt) ~ educ * light_30_log + age + female + 
                     urban + unempl + ethfair + media +
                     conflict + ea_service + ea_infrastructure + as.factor(cyear) +
                     (1 | EA), data = petro, family=cumulative("logit"),
                   chains=4, seed=123, iter=2000, warmup=1000)
#saveRDS(liv30_petro, "fit/liv30_petro.Rds")
mytable(liv30_petro, livingcond = T)

## Varying sizes of local economic environments---- 
# Dynamic Radii
ecodyn <- brm(probecon | weights(combinwt) ~ educ * light_dyn_log + age + female + 
                urban + unempl + ethfair + media +
                conflict + ea_service + ea_infrastructure + as.factor(cyear) +
                (1 | EA), data = dat, family=binomial("logit"),
              chains=4, seed=123, iter=2000, warmup=1000)
#saveRDS(ecodyn, "fit/econdyn.Rds")
mytable(ecodyn, livingcond = F)

livdyn <- brm(as.ordered(plivcond)  | weights(combinwt) ~ educ * light_dyn_log + age + female + 
                urban + unempl + ethfair + media + 
                conflict + ea_service + ea_infrastructure + as.factor(cyear) +
                (1 | EA), data = dat, family=cumulative("logit"),
              chains=4, seed=123, iter=2000, warmup=1000)
#saveRDS(livdyn, "fit/livdyn.Rds")
mytable(livdyn, livingcond = T)

# 10km radius
econ10 <- brm(probecon | weights(combinwt) ~ educ * light_10_log + age + female + 
                urban + unempl + ethfair + media + 
                conflict + ea_service + ea_infrastructure + as.factor(cyear) +
                (1 | EA), data = dat, family=binomial("logit"),
              chains=4, seed=123, iter=2000, warmup=1000)
#saveRDS(econ10, "fit/econ10.Rds")
mytable(econ10, livingcond = F)

liv10 <- brm(as.ordered(plivcond)  | weights(combinwt) ~ educ * light_10_log + age + female + 
               urban + unempl + ethfair + media + 
               conflict + ea_service + ea_infrastructure + as.factor(cyear) +
               (1 | EA), data = dat, family=cumulative("logit"),
             chains=4, seed=123, iter=2000, warmup=1000)
#saveRDS(liv10, "fit/liv10.Rds")
mytable(liv10, livingcond = T)

## 50km radius
econ50 <- brm(probecon | weights(combinwt) ~ educ * light_50_log + age + female + 
                urban + unempl + ethfair + media + 
                conflict + ea_service + ea_infrastructure + as.factor(cyear) +
                (1 | EA), data = dat, family=binomial("logit"),
              chains=4, seed=123, iter=2000, warmup=1000)
#saveRDS(econ50, "fit/econ50.Rds")
mytable(econ50, livingcond = F)

liv50 <- brm(as.ordered(plivcond)  | weights(combinwt) ~ educ * light_50_log + age + female + 
               urban + unempl + ethfair + media + 
               conflict + ea_service + ea_infrastructure + as.factor(cyear) +  (1 | EA),
             chains=4, seed=123, iter=2000, warmup=1000,
data = dat, family=cumulative("logit"))
#ßsaveRDS(liv50, "fit/liv50.Rds")
mytable(liv50, livingcond = T)

## Non-hierarchical, frequentist models----
r1 <- felm(probecon ~ educ*light_30_log + age + female + urban + unempl + ethfair + media + ea_service + ea_infrastructure + conflict | raster50 + cyear | 0 | EA, data = dat)
m1 <- felm(probecon ~ educ*light_30_log + age + female + urban + unempl + ethfair + media | EA + cyear | 0 | EA, data = dat)

r2 <- felm(plivcond ~ educ*light_30_log + age + female + urban + unempl + ethfair + media +ea_service + ea_infrastructure + conflict | raster50 + cyear | 0 | EA, data = dat)
m2 <- felm(plivcond ~ educ*light_30_log + age + female + urban + unempl + ethfair + media | EA + cyear | 0 | EA, data = dat)

# Table 18
stargazer::stargazer(r1, m1, r2, m2, 
                     order = c(1,2,12,3,4,5,6,7,8,9,10,11),
                     dep.var.labels.include = TRUE,
                     dep.var.labels = c("Labor market risks","Living conditions"),
                     keep.stat = c("n","rsq"),
                     covariate.labels=c("Education", "Local development", "Education * Local development", 
                                        "Age", "Female", "Urban", "Unemployed", "Ethnic grievance", "Media consumption",
                                        "Service (EA)", "Infrastructure (EA)", "Conflict (EA)"),
                     add.lines = list(c("Country-Year FE", "Yes", "Yes", "Yes", "Yes"),
                                      c("Raster FE", "Yes", "", "Yes", ""),
                                      c("EA FE", "","Yes", "", "Yes")),
                     star.cutoffs = c(.1, .05), 
                     star.char = c("\\dagger", "*"),
                     notes = "Standard errors clustered on survey cluster; $^{\\dagger} p <$ 0.1, $^*p <$ 0.05",
                     notes.append = F,
                     table.placement = "ht!", column.sep.width = "1pt",
                     digits = 2,
                     title = "",
                     label = "tab:linear",
                     type="text")

